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I.  INTRODUCTION 


Systems  are  modelled  In  order  to  understand  and  explain  them  better 
and  as  a  prelude  to  action.  Aircraft  dynamics*  for  example*  may  be 
Identified  so  that  better  designs  can  be  made,  or  so  that  adaptive  control 
actions  can  be  taken.  Our  attention  in  this  study  has  been  directed  at 
aspects  of  estimation  and  Identification  that  are  connected  with  system 
understanding. 

This  study  has  been  aimed  at  three  problem  areas:  multistage 
modeling*  estimation,  and  Identification  algorithms  for  dynamical  systems: 
parametric  Identification  algorithms  for  estimation  of  Input-output 
response;  and,  state  estimation  and  parameter  Identification  for  a  new 
class  of  models  —  causal  functional  equations  ~  which  describe  wave 
propagation  In  layered  media  systems. 

The  first  problem  Is  concerned  with  developing  estimation  algorithms 
both  for  parameter  and  state  estimation  that  are  recursive  In  the  dimen¬ 
sion  of  the  parameter  vector  or  state  vector.  Such  algorithms  will  find 
utility  in  system  modeling  work,  where  model  dimension  Is  often  a  variable. 

The  second  problem  Is  concerned  with  Identifying  a  linear  system's 
Input-output  response  rather  than  a  parametric  representation  for  the 
sytem.  The  basic  Idea  is  that  parameter  estimation  techniques  often  give 
very  poor  estimates  of  a  system's  response;  hence,  we  are  Interested  In 
new  algorithms  that  lead  to  good  estimates  of  system  Input-output  response. 

The  third  problem  Is  concerned  with  developing  whole  new  theories  of 
state  estimation  and  parameter  Identification  for  causal  functional  equa¬ 
tions.  These  equations  are  continuous-time,  linear,  time-invariant  and 
contain  multiple  time  delays.  They  do  not  contain  derivatives  or  Integrals, 
and  no  literature  apparently  exists  for  them.  Causal  functional  equations 
are  applicable  to  diverse  areas  such  as  reflection  seismology,  transmission 
lines,  speech  processing,  optical  thin  coatings  and  EM  problems. 


II.  SUMMARY  OF  RESEARCH  EFFORT 


We  sunmarlze  our  research  activities  for  the  three  problem  areas 
mentioned  In  Section  I  In  Paragraphs  A.  B,  and  C  below.  Some  tangential 
results  are  also  summarized  In  Paragraph  0. 

A.  Multistage  Algorithms 

1.  Multistage  least-squares  parameter  estimation  algorithms  have 
been  extended  to  multistage  weighted  least-squares  algorithms  [1].  This 
extension  permits  one  to  do  multistage  unbiased  minimum  variance  estima¬ 
tion  of  parameters  In  a  linear  model,  which  broadens  the  multistage 
philosophy  to  more  than  one  estimation  technique. 

2.  Multistage  least-squares  algorithms  have  been  applied  to  a 
number  of  Interesting  problems  associated  with  Identification  of  a  sampled 
Impulse  response  [1].*  Suppose  a  test  signal,  u(t).  Is  applied  at  t»tQ 

to  a  linear,  time-invariant,  causal,  but  unknown  system  whose  output, 
y(t),  and  Input  are  measured.  The  unknown  Impulse  response  is  to  be 
Identified  using  sampled  values  of  u(t)  and  y(t).  For  such  a  system 

y(tk)  »/  k  w(-c)u(tk  -  T)dx  ;  (1) 

t0 

or.  If  we  assume  that  (a)  w(t)  »  0  for  all  t  2  tw,  (b)  [tQ,  t^]  Is 
divided  Into  N  equal  Intervals,  each  of  width  AT,  and  (c)  for 
T€CtM.  tj],  w(t)  -w(ti_1),  and  u(t  -  t)  -  u(t  -  t^),  then 


♦Publications  under  this  grant  are  listed  in  Section  III.  References 
[e.g.[A]]  are  listed  at  the  end  of  the  report. 


(2) 


y(V  *■  ^  wl(ti-l)u(tk  "  *1-1* 


where 


*1^*1-1*  *  AT  w(ti_1)  .  (3) 

It  Is  straightforward  to  Identify  the  N  unknown  parameters  w^tg), 
w-|(t, *|(tN_i)  via  least-squares;  however,  for  N  to  be  known  t 
must  be  accurately  known  and  AT  must  be  chosen  judiciously. 

In  actual  practice  t^  Is  not  known  that  accurately  so  that  N  may  have 
to  be  varied.  Multistage  least-squares  estimators  (LSE's)  have  been  used 
to  handle  this  situation  In  a  computationally  efficient  manner. 

Sometimes,  AT  can  be  too  coarse  for  certain  regions  of  time.  In  which 
case  significant  features  of  w(t),  such  as  a  ripple,  may  be  obscured.  In 
this  situation,  we  would  like  to  Hzoom-1nH  on  those  Intervals  of  time  and 
re-dlscretlze  y(t)  just  over  those  Intervals,  thereby  adding  more  terms  to 
yU^).  Me  have  shown  how  to  use  multistage  LSE's  to  zoom-in  on  regions  of 
time  In  a  very  effective  manner. 

3.  Computational  requirements  for  multistage  LSE's  have  been  studied 
[2].  Our  approach  was  to  count  the  total  number  of  multiplications  asso¬ 
ciated  with  our  different  multistage  algorithms.  Total  number  of  multipli¬ 
cations  Is  a  good  Indicator  of  total  computation  time.  We  have  compared 
the  computation  time  requirements  for  our  multistage  algorithms  with  one- 
stage  (conventional)  LSE  algorithms  from  two  different  points  of  view. 

In  order  to  Illustrate  the  different  points  of  view,  recall  that  by 
means  of  n-stage  LSE's  It  Is  possible  to  simultaneously  obtain  least- 
squares  estimates  for  parameters  In  n  models.  According  to  the  first  point 
of  view  we  treat  the  n-stage  LSE  as  a  recursive  procedure  for  obtaining 
just  the  least-squares  estimate  of  n-vector  e;  hence,  the  computation  time 
required  to  obtain  e  by  means  of  the  n-stage  LSE  should  be  compared  with 

“  A 

the  computation  time  required  to  obtain  e  by  means  of  a  one-stage  LSE.  By 
this  first  point  of  view  we  are  comparing  computation  time  requirements  on 
an  absolute  basis.  According  to  the  second  point  of  view  we  treat  the 
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n- stage  LSE  as  a  recursive  procedure  for  obtaining  least-squares  estimates 
for  n  different  linear  models;  hence,  the  computation  time  required  to 
obtain  9  by  means  of  the  n-stage  LSE  should  be  compared  with  the  computa¬ 
tion  time  required  to  obtain  least-squares  estimates  for  each  one  of  the 
n  linear  models,  each  estimate  being  obtained  by  a  one-stage  LSE.  By  the 
second  point  of  view  we  are  comparing  computation  time  requirements  on  the 
basis  of  equal  Information. 

The  following  conclusions  have  been  reached  (quantitative  conclusions 
are  given  In  [2]):  For  batch  data  processing.  It  Is  always  more  efficient 
to  use  multistage  LSE's.  This  Is  true  both  on  an  absolute  and  equal  Infor 
mat ion  basis.  For  sequential  data  processing,  conclusions  depend  on 
whether  we  choose  to  use  a  "standard"  LSE  or  a  "stabilized"  LSE.  Often, 
for  purposes  of  accuracy,  one  must  use  a  stabilized  LSE,  In  which  case  It 
Is  again  always  more  efficient  to  use  multistage  LSE's.  If,  on  the  other 
hand,  accuracy  Is  not  so  Important,  then  It  Is  usually  more  efficient  to 
use  one-stage  standard  LSE's  than  multistage  standard  LSE's. 

4.  Frledland's  [A]  bias  filtering  technique,  which  Is  an  example  of 
a  multistage  Kalman/Bucy  (K/B)  filter  for  a  very  special  plant  (or  transi¬ 
tion)  matrix  and  process  covariance  matrix,  has  been  extended  to  the 
following  class  of  nonlinear  systems  [3]: 


x(t)  *  f[x(t),  t]  +  A-jb  +  w(t) 

(4a) 

y(t)  «  h[x(t)»  t]  +  A2b  +  v(t). 

(4b) 

In  this  system,  states  x  enter  nonllnearly;  but,  biases  enter  linearly. 
Signals  w(t)  and  v(t)  are  zero  mean  white  noise  processes. 

5.  In  our  efforts  to  generalize  Frledland's  [A]  decomposition  to  the 
situation  of  a  bias  vector  described  by  a  first-order  Markov  process,  we 
have  had  to  rederive  his  decomposition  by  a  route  different  from  the  one 
taken  by  him.  His  approach  was  to  augment  the  bias  states  to  the  dynamical 
x  states  and  to  show  that  the  Rlccati  equation  for  the  augmented  system 
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could  be  decomposed  so  that  x.(k)  can  be  computed  via  the  following 
equation: 

x(k)  -  x(k)  +  Vx(k)b(k)  (5) 

where  x(k)  Is  the  bias-free  estimate  of  x(k),  computed  as  If  no  biases 

“  A  " 

were  present,  b(k)  Is  the  optimum  estimate  of  the  bias,  and  Vw(k)  Is  a 
matrix,  that  Is  computed  recursively,  which  blends  the  estimates  of  x(k) 

A  A 

and  b(k)  together  to  give  x(k),  the  bias  corrected  estimate  of  x(k). 

Our  approach  [4,5],  which  Is  constructive  In  nature.  Is  to  assume  the 
existence  of  the  decomposition  In  Eq.  (5),  for  which  we  require  that 

A  A 

x(k)  and  b(k)  be  unbiased  estimators  of  x(k)  and  b(k),  and,  that  the 
trace  of  the  estimation  error  covariance  matrix  for  x(k)  be  minimum. 

A  " 

We  then  demonstrate  that  the  resulting  x(k)  Is  Indeed  the  minimum 
variance  estimator  of  x(k);  hence,  our  assumed  decomposition  Is  valid. 

Besides  obtaining  all  of  Frledland's  results  (and,  giving  some  new 
physical  meaning  to  same  of  the  matrix  quantities  such  as  Vv(k)),  we 
also  show  that  b(k)  In  Eq.  (5)  Is  a  minimum  variance  estimator  of  b  for 
an  auxiliary  parameter  estimation  problem;  hence,  we  have  demonstrated 
that  It  Is  possible  to  add  or  remove  elements  from  b  so  that  x(k)  and 
the  modified  bias  vector  can  be  reestimated  without  having  to  redo  all 
the  calculations.  This  Is  Important  In  applications  where  the  bias 
states  are  used  to  model  constant  Instrumentation  error  sources  [B];  for. 
It  can  happen  that  not  all  of  the  error  sources  are  significant,  so  that 
some  of  them  should  be  deleted  from  the  final  filter,  or,  that  signifi¬ 
cant  error  sources  may  have  been  Initially  neglected,  and  should  be 
Included  In  the  final  filter. 

6.  Washburn  [6,7]  has  generalized  Frledland's  [A]  bias  estimation 
technique  to  partitioned  dynamical  systems.  In  the  general  case,  the 
calculations  are  of  the  dimension  of  the  overall  system,  so  that,  except 
for  some  special  but  Important  cases,  there  are  no  computational  advan¬ 
tages  to  the  multistage  approach.  Those  special  cases,  where  there  does 
appear  to  be  computational  advantages  for  the  multistage  approach  are: 
colored  noise  and  weak  coupling  between  the  partitioned  systems. 
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The  general  results  are  Important  In  themselves,  since  they  provide 
the  theory  for  a  particular  decomposition  of  the  optimal  state  estimator 
for  a  system  of  possibly  large  dimension  (l.e.,  a  large  scale  system). 

This  decomposition  gives  added  insight  Into  the  structure  and  performance 
of  the  minimum  variance  unbiased  estimator.  In  addition,  the  methodology 
of  proof  for  this  multistage  decomposition  provides  a  means  for  Investi¬ 
gating  other  decompositions  of  Interest. 

B.  Identification  From  Input/Output  Data 

1.  Several  existing  Identification  algorithms  have  been  analyzed 
to  determine  how  well  the  resulting  models  match  the  Input/output 
response  of  the  true  system,  especially  with  respect  to  stability.  It 
has  been  found  that  In  the  Identification  of  autoregressive  processes 
the  standard  least-squares  algorithm,  (which  Is  also  maximum  likelihood) 
gives  good  results.  This  would  be  expected  from  the  properties  of 
maximum  likelihood  methods,  and  the  Inherent  numerical  robustness  of  the 
algorithm.  This  robustness  Is  due  to  the  fact  that  a  Toeplltz  matrix  Is 
Inverted  which  Is  generally  positive  definite  and  well -conditioned. 

However,  In  the  Identification  of  autoregressive  moving  average  processes 
a  number  of  algorithms  first  estimate  the  Markov  parameters  (e.g.,  [C], 

[0])  and  then  Invert  the  associated  Hankel  matrix.  It  has  been  found 
that  these  algorithms  do  not  perform  well  In  their  reproduction  of  the 
input/output  response  for  two  reasons.  Firstly,  because  only  a  small 
number  of  Markov  parameters  Is  used,  there  Is  high  sensitivity  to  errors, 
and  secondly,  the  Hankel  matrix  may  become  almost  singular  (Glover  [E]). 
These  results  have  been  verified  by  simulation  for  the  method  given  In  [C]. 

2.  Preliminary  work  on  Input  design  for  identification  [8]  has 
Indicated  that  the  choice  of  design  criterion  can  be  critical.  An 
illustrative  example  was  produced  where  the  trace  of  the  Information 
matrix  was  maximized  subject  to  an  Input  energy  constraint  as  follows. 

There  are  two  Identifiable  parameters  to  be  identified  and  the  available 
Input  energy  must  be  shared  between  both.  The  above  criterion  however 
indicated  that  all  the  Input  energy  should  be  channelled  Into  Identifying 
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the  more  easily  Identified  parameter.  This  resulted  In  one  parameter 
being-  Identified  very  well  and  no  Information  being  obtained  on  the  other 
parameter,  which  Is  clearly  not  desired. 

3.  Several  deterministic  Identification  problems  have  been  solved 
[9].  The  general  problem  Is:  given  a  finite  set  of  finite  length  Input/ 
output  sequences  from  an  unknown  discrete-time  system,  derive  a  minimal 
order  linear  dynamical  model  with  appropriate  Initial  conditions  that 
would  produce  the  observed  Input-output  sequences.  Previous  approaches 
to  this  problem  (e.g.  [F])  have  been  based  on  Involved  matrix  manipula¬ 
tions  and  have  required  certain  restrictive  assumptions  on  the  data.  The 
new  method  Is  based  on  the  geometric  theory  of  linear  systems  and  Is 
completely  general.  In  that  an  arbitrary  number  of  Input/output  sequences 
of  arbitrary  length  are  allowed,  and  the  Initial  conditions  can  be  either 
zero  or  free.  The  resulting  algorithm  Is  simply  stated  In  terms  of 
recursive  subspace  equations  and  Includes  tests  for  the  Identlf lability 
of  the  data.  It  has  further  been  shown  that  several  other  problems  In 
linear  systems  (e.g.  observers  and  exact  model  matching)  can  be  formulated 
within  the  same  framework. 

4.  Substantial  progress  has  been  made  In  Implementing  and  evaluating 
approximate  realization  methods  (l.e.,  the  approximate  realization  of 
given  Impulse  response  sequence  by  a  low  order  linear  system).  Given 
exact  data  there  are  a  large  number  of  different  algorithms  that  will 
produce  exact  realizations,  and  given  approximate  data  the  same  methods 
can  be  approximately  Implemented  to  produce  approximate  realizations. 

This  Is  similar  In  spirit  to  an  equation  error  method  and  a  selection  of 
six  approximate  realization  methods  of  this  type  have  been  Implemented 
and  have  been  evaluated  using  experimental  Impulse  response  data. 

C.  Modeling,  Estimation,  and  Identification  of  Layered  Media  Systems 

1.  We  have  developed  time-domain  state  space  models  for  lossless 
layered  media  which  are  described  by  the  wave  equation  and  boundary 
conditions  [10-13].  Our  models  are  for  non-equal  one-way  travel  times; 
hence,  they  are  more  general  than  existing  models  of  layered  media  which 
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are  usually  for  layers  of  equal  one-way  travel  times.  Full  state  models* 
which- Involve  2K  states  for  a  K-layer  media  systems*  as  well  as  half¬ 
state  models*  which  Involve  only  K  states  have  been  developed  and  related. 
Certain  transfer  functions*  which  appear  In  the  geophysics  literature  In 
connection  with  models  of  layered  media  with  equal  travel  times  have  been 
generalized  to  the  situation  of  non-equal  travel  times.  Our  state  space 
models  represent  a  new  class  of  equations,  causal  functional  equations, 
for  which  we  have  not  been  able  to  find  any  literature.  These  equations 
are  continuous-time,  linear,  and  contain  multiple  time-delays.  Their 
Impulse  response  Is  an  Infinite  sequence  of  non-unlformly  spaced  Impulse 
functions. 

2.  We  have  proven  the  truth  [11,14]  of  the  following  decomposition 
of  the  solutions  to  the  lossless  wave  equation  In  layered  media:  the 
complete  output  from  a  K-layer  media  system,  which  Is  comprised  of  the 
superposition  of  primaries,  secondaries,  tertlarles,  etc.,  can  be  obtained 
from  a  single  state  space  model  of  order  2K  —  the  complete  model  —  or 
from  an  Infinite  number  of  models,  each  of  order  2K,  the  output  of  the 
first  of  which  Is  just  the  primaries,  the  output  of  the  second  of  which 
Is  just  the  secondaries,  etc.  This  decomposition  of  the  solution  to  the 
lossless  wave  equation  Into  physically  meaningful  constituents  (l.e., 
primaries,  secondaries,  etc.)  Is  called  a  canonical  Bremmer  Series  de¬ 
composition,  after  Bremmer,  who  in  1951  established  a  similar  decomposi¬ 
tion  [6]. 

In  many  geophysical  situations,  where  reflection  coefficients  are 
quite  small,  the  decomposition  can  be  truncated  after  secondaries  or 
tertlarles;  hence.  It  also  represents  a  way  to  approximate  the  solution 
to  the  wave  equation. 

We  have  made  connections  [15]  between  our  state  space  models  and  the 
Integral  equations  given  by  Bremmer  [G]  for  generating  the  partial 
residuals. 

We  have  shown  how  to  go  from  Brenner's  Integral  equations  to  our 
state  equations  by  assuming  a  medium  with  a  wave  number  that  has  finite 
jumps  (discontinuities)  which  occur  at  the  Interfaces  and  that  Is  constant 
within  a  layer.  These  assumptions  for  wave  number  are  associated  with 
what  we  mean  by  a  horizontally  layered  homogeneous  earth. 


-8- 


We  have  also  demonstrated  that  Brenner's  Integral  equations  can  be 
obtained  by  the  W.  K.  B.  method  which  gives  approximate  solutions  to 
second-order  differential  equations  [H].  This  justifies  earlier  claims 
[11*14]  that  the  Brenner  series  decomposition  can  be  used  to  approximate 
the  complete  solution  to  the  lossless  wave  equation  by  truncating  that 
decomposition  after  a  small  number  of  terms. 

3.  We  have  developed  [16]  a  general  theory  for  describing  reinforced 
events  between  multiple  reflections  in  lossless  layered  media,  which  are 
described  by  the  wave  equation  and  boundary  conditions  (e.g.*  horizontally 
stratified  nonabsorptive  earth  with  vertically  traveling  plane  congres¬ 
sional  waves). 

Reinforcements  occur  whenever  two  or  more  multiple  reflections  from 
different  paths  inside  the  media  arrive  at  the  surface  at  the  same  time 
so  that  they  add  (positively  or  negatively)  together.  Those  reinforce¬ 
ments  occur  regardless  of  what  the  travel  time  is  in  each  layer,  and 
distort  the  appearance  of  a  seismogram;  for,  they  lead  one  to  believe  that 
a  significant  event  has  occurred  by  the  appearance  of  a  large  amplitude 
segment  of  the  seismogram,  whereas,  in  reality,  that  large  event  is  a  sum 
of  (many)  smaller  events. 

Our  general  theory  Is  applicable  to  a  K-layer  media  system  with  non- 
uniform  travel  times  and  gives  Information  about  the  exact  location  in 
time,  number,  and  amplitude  of  reinforced  events  for  n-arles  (i.e., 
secondaries,  tertlaries,  etc.),  where  n  *  1,  2,  3  ...  .  The  starting 
point  for  the  development  of  this  theory  is  Mendel's  Bremmer  series 
decomposition  [11,14]  and  the  operator  description  of  state  space  models 
of  layered  media  [12]  by  means  of  which  n-ary  reflections  (where  n  ■  1, 

2,  3,  ...)  are  generated  and  analyzed  separately  and  related  to  each 
other.  The  two  most  significant  multiple  reflections,  secondaries,  and 
tertlaries,  have  been  studied  extensively.  We  have  demonstrated  that 
not  only  do  reinforcements  occur  between  the  same  kind  of  multiple 
reflections  (e.g.,  between  secondaries),  but  that  reinforcements  also 
occur  across  different  kinds  of  multiple  reflections  (e.g.,  between 
secondaries  and  tertlaries). 


4.  Because  our  causal  functional  state  space  models  for  a  layered 
media  system  represent  a  new  class  of  equations,  we  have  had  to  study  the 
computer  simulation  of  these  equations.  Two  computational  methods  have 
been  considered  [17].  In  the  first  approach,  we  discretized  the  time 
axis  and  Inserted  states  of  Intermediate  delays,  to  arrive  at  a  set  of 
standard  finite-difference  equations.  For  our  particular  system,  matrix 
multiplications  can  be  reduced  to  simple  scalar  multiplications.  In  the 
second  approach,  we  defined  mapping  rules  for  the  transformation  of  states 
at  an  Interface,  and  kept  a  state  reference  table  for  look-up  and  branching. 
The  procedure  Is  similar  to  ray-tracing.  Several  experiments  have  been 
performed  to  show  the  trade-off  between  storage  requirement  and  CPU  time- 
spent  for  the  two  methods. 

5.  We  have  developed  a  procedure  for  extracting  reflection  coeffi¬ 
cients  from  noisy  data  [18]  which  we  feel  Is  a  substantial  generalization 
of  similar  procedures  which  have  been  reported  In  the  literature  ([I]  for 
example).  Associated  with  these  earlier  procedures  are  Standard  Assump¬ 
tions  and  Steps  which  Include  requirements  that  the  data  be  noise  free 
and  that  the  observed  seismic  data  be  deconvolved.  Our  procedure  avoids 
these  restrictive  requirements.  Furthermore,  our  procedure  totally  avoids 
the  concepts  of  z-transforms,  minimum  phase,  spectral  factorization,  etc., 
which  appear  in  the  literature  on  this  subject. 

6.  An  Important  special  case  of  a  causal  functional  equation  (CFE), 

occurs  when  all  one-way  travel  times  are  equal.  In  this  case  the  uniform 
CFE  (UCFE)  is  ^ 

x(t  +  t)  »  A  x_(t)  +  m(t)  (6) 


with  Initial  values 

x(o)  o  €  [0,  t)  £  T  (7) 

Recognizing  that  any  time  t  (t  €  R)  can  be  expressed  as 

t  *  t'  +  Mr  where  t'6  J  and  M  is  an  Integer  (8) 
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we  have  shOMn  [19]  that  the  solution  to  UCFE  (6)  Is 


x  [t'  +  (k+1)T]  -  Ak+1  x(f)  +  S  A*"1  m(t*  +  It)  (9) 

1*0 

where  t'  6  T*  k  *  0,  1,  2,  ...,  and  t  ■  t'  +  (k+l)t. 

Equation  (9)  explicitly  shows  how  the  state  at  any  time  t  *  t'  +  (k+l)r 
depends  on  an  Initial  condition  x_(t')  and  the  Input  m.  It  Is  of  Interest  to 
note  that  x(t)  depends  only  on  a  single  element  of  the  Initial  values  x(o) 

(o  €  j)»  namely  x(t'),  and  a  finite  number  of  point  values  of  m.  This  shows 
that  the  solution  to  the  uniform  causal  functional  state  equation,  although 
continuous- time  In  nature,  derives  Its  values  In  a  discrete-time  fashion  for 
a  given  fixed  value  of  t*  €  J.  Of  course,  there  are  an  uncountable  number 
of  points  In  T;  hence  we  can  imagine  x(t)  as  being  generated  by  an  uncountable 
number  of  discrete- time  systems. 

When  we  simulate  our  results  on  a  digital  computer,  computations  are 
made  every  T  sec.  at  discrete  time  points.  Consequently,  on  a  digital  com¬ 
puter,  x(t)  Is  generated  by  a  finite  number  of  discrete-time  systems  which 
operate  In  parallel. 

7.  We  have  derived  [19]  the  minimum-variance  state  estimator  for  UCFE 
x(t  +  t)  *  A  x(t)  +  B  m(t)  +  w(t)  (10) 

and  Its  associated  measurement  equation 

X(t)  *  H  x(t)  +  n(t)  (11) 

a 

Let  x(t)  denote  the  minimum-variance  estimate  of  x(t)  based  on  all 
measurements  In  {^(X):  0  <  X  s  t,  t  €  R).  We  have  shown  that  for  any 
fixed  t'  6  Tm  [0,  t),  x(t),  where  t  *  t'  +  Mr  (M  *  1,  2,  ...),  Is  given 
by  the  usual  discrete-time  Kalman  filter  equations  (Ref.  J,  for  example) 

A 

with  t1  considered  the  Initial  starting  time.  Of  course,  to  obtain  x(t) 
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for  all  t  6  R  we  would  need  an  uncountable  number  of  discrete-time  Kalman 
filters;  but.  Imposing  a  mesh  on  j  (with  grid  size  equal  to  data  sampling 
rate)  leads  to  a  finite  number  of  Kalman  filters  which  operate  In  parallel. 
To  the  best  knowledge  of  the  author  this  Is  the  first  estimation  theory 
result  that  has  led  to  a  natural  form  of  parallel  data  processing. 

8.  We  have  developed  an  extended  minimum-variance  estimator  for 
simultaneous  estimation  of  states  and  parameters  (l.e.,  reflection  co¬ 
efficients)  In  a  UCFE  [20].  Simulation  results  were  very  disappointing. 
Reasons  for  the  disappointing  results  are  explained  in  [20]. 

9.  A  simple  Inverse  filter  has  been  developed  [21]  to  suppress 
multiple  reflections  from  a  normal  Incidence  synthetic  seismogram.  The 
filter  was  developed  by  means  of  Mendel's  Bremmer  Series  Decomposition 
[14]  and  the  operator  description  of  state  space  model  of  layered  media 
[12],  and  Is  given  In  terms  of  z-transforms  as 

Y-j  (z)  «  - 
1 


In  this  equation  Y(z)  Is  the  synthetic  seismogram  measured  by  a  sensor 
located  on  the  surface.  That  sensor  receives  reflected  signals  from  a 
layered  media  due  to  the  normal  Incident  Input  M(z)  applied  at  the  same 
surface.  Parameter  rQ  Is  the  reflection  coefficient  of  the  surface. 

Y,(z)  Is  the  output  of  the  filter;  It  consists  of  the  primary  reflection 
portion  of  the  seismogram  and  some  residual  terms;  l.e.,  Y^z)  *  Y^(z)+y(z)* 
In  general,  the  residual  terms  y{z)  Is  relatively  quite  small,  and  Y^z) 

Is  a  good  approximation  of  Y^z).  This  filter  Is  especially  effective 
when  rQ  Is  relatively  large  (as  In  most  geophysical  situations)  In  which 
case  y(z)  Is  almost  negligible  compared  with  Yj(z).  The  filter  requires 
knowledge  of  the  Input  waveform  M(z),  surface  reflection  coefficient,  Tq, 
and  measured  seismogram  Y(z).  Observe  that  (12)  represents  a  nonlinear 
processing  of  the  seismogram  Y(z). 


10.  We  have  extended  our  normal  Incidence  state  space  model  to  the 
non-normal  Incidence  case  [22].  The  non-normal  Incidence  (NNI)  state 
space  model  Is  structurally  the  same  as  the  normal  Incidence  state  space 
model  except  that  It  has  twice  as  many  state  variables.  Because  of  mode 
conversion  In  non-normal  Incidence,  the  scalar  upgolng  and  downgoing  waves 
and  travel  times  In  each  layer  as  well  as  reflection  and  transmission  co¬ 
efficients  In  each  Interface  are  replaced  by  a  vector  of  upgolng  and  down¬ 
going  waves,  a  vector  of  travel  times,  and  matrices  of  reflection  and 
transmission  coefficients. 

With  this  NNI  model,  we  are  able  to  generate  synthetic  seismograms 
for  a  plane  wave  source,  and  more  Importantly,  for  a  two-dimensional 
point  source. 

11.  We  have  developed  a  maximum-likelihood  procedure  for  estimating 
both  the  reflection  coefficients  and  one-way  travel  times  [23]  for  a 
lossless  layered  media  system  in  which  the  layers  are  non-equal ly  spaced 
(In  time).  It  uses  a  state  space  model  as  Its  starting  point,  one  that 
Is  more  general  than  (6)  since  now  t  Is  different  for  each  layer  (l.e., 

x1  t  t2  f  x3  f  ...  Tk  t  t).  The  only  source  of  uncertainty  Is  measurement 
noise,  n(t).  Maximum-likelihood  estimates  of  the  parameters  are  obtained 
In  a  layer-recursive  format.  In  essence,  first  r1  and  t-j  are  determined 
and  layer  1  Is  stripped  away;  then  r2  and  t2  are  determined  and  layer  2 
Is  stripped  away;  etc.  To  the  best  of  our  knowledge,  this  Is  the  first 
time  that  both  reflection  coefficients  and  travel  times  have  been  simul¬ 
taneously  estimated  In  an  optimal  manner. 

12.  We  have  surveyed  approaches  to  solving  inverse  problems  for 
lossless  layered  media  systems  [24]. 

13.  A  Kunetz  equation  Is  often  used  as  the  starting  point  In  the 
development  of  solutions  for  the  Inversion  of  one-dimensional,  noise-free, 
normal  Incident  seismograms,  for  which  |rg|  *  1  (rg  Is  the  surface 
reflection  coefficient).  We  have  demonstrated  a  need  for  a  Kunetz-type 
equation  In  which  filtered  signals  can  be  used,  so  that  noise  effects 
(which  are  always  present  In  real  data)  can  be  reduced.  Furthermore,  we 


have  shown  that  an  Infinite  numberof  Kunetz-type  equations  exist  for  the 
lossless  wave  equation  In  layered  media.  Finally,  we  have  shown  that  It 
Is  Indeed  valid  to  formulate  and  solve  the  Inverse  problem  using  filtered 
signals  [25]. 

14.  We  have  developed  a  unified  procedure  for  estimating  reflection 
coefficients  of  lossless  layered  media  systems  from  noise  seismic  data 
[26].  Our  procedure  Is  a  generalization  of  work  described  In  [18]  to  the 
cases  of  source  and  sensor  either  at  the  surface  or  In  the  first  layer 
(e.g.,  a  water  layer).  We  have  handled  these  cases  simultaneously,  since, 
conceptually  at  least,  they  are  similar. 

Our  procedure  Is  to:  (1)  write  state  equations  which  describe  the 
temporal  evolution  of  all  upgolng  and  downgoing  signals;  (2)  derive  an 
autoregressive-like  equation  for  basement  signal  dK+1(t);  (3)  establish  a 
parameter  estimation  problem  for  estimating  the  coefficients  In  the 
d|(+l(t)  e(Juat1°n’  and  generate  the  Normal  equations  which  result  from 
solution  of  that  problem;  (4)  express  the  Normal  equations  In  terms  of 
measurable  signals,  through  use  of  a  Kunetz  equation;  and  (5)  study 
conditions  for  which  the  Normal  equation  can  be  solved,  as  well  as  proper¬ 
ties  of  the  solution. 

Our  procedure  Is  constructive  and  during  the  study  of  Its  five  steps 
we  have  provided  answers  to  the  following  questions:  (1)  Why  do  we  direct 
our  attention  at  an  equation  for  the  unmeasurable  basement  signal,  dK+1(t), 
rather  than  the  measurable  seismogram  signal?;  (2)  What  Is  the  real  role 
of  the  Kunetz  equation?;  (3)  Where  In  the  procedure  Is  It  necessary  to 
distinguish  the  cases  of  source  and  sensor  at  the  surface  versus  source 
and  sensor  in  the  first  layer?;  (4)  Why  Is  deconvolution  necessary?;  (5) 

Why  Is  the  solution  limited  to  cases  for  which  the  amplitude  of  the  surface 
reflection  coefficient  Is  unity?;  and  (6)  Why  is  the  solution  so  sensitive 
to  noise? 

15.  We  have  sunmarlzed  our  work  on  uniform  causal  functional 
equations  (UCFE)  [27],  (10)  and  (11).  Because  a  UCFE  Is  Isomorphic  to  an 
uncountable  number  of  discrete-time  systems,  each  one  initialized  at 

t'  6  T  *  [0,  t],  we  have  been  led  to  a  very  simple  proof  of  a  widely  used 
fact  that  UCFE  (10)  Is  properly  Initialized  by  x(o),  a  €  J". 


16.  We  have  continued  work  on  Hablbl-Ashrafl's  [23]  maximum-likelihood 
procedure  for  estimating  reflection  coefficients  and  detecting  one-way 
travel  times  [28]  (In  all  of  our  other  studies  travel  time  t  1$  assumed 
known,  and  all  layers  have  the  same  travel  time).  We  are  especially 
Interested  In  making  connections  between  Ms  results  and  our  suboptlmal 
results  [26],  when  Ms  results  are  specialized  to  the  case  of  uniform  travel 
time. 


17.  We  have  developed  a  Brenner  series  decomposition  for  the  NNI  case 
and  have  also  designed  approximate  filters  for  suppression  of  multiple 
reflections  [29,30].  These  filters  are  based  on  the  algebraic  structure  of 
the  NNI  Brenner  series  decomposition. 

D.  Miscellaneous 


1.  It  Is  common  In  system  modeling  that  a  number  of  parameters  are 
not  known  precisely  but  will  be  determined  later  from  empirical  data  or  by 
subsequent  design  decisions.  It  Is  therefore  Important  to  study  such 
structured  systems,  and  In  [31]  we  have  derived  necessary  and  sufficient 
conditions  for  the  structural  controllability  of  multi -Input  structured 
linear  systems.  The  methods  used  In  [31]  are  substantially  simpler  than 
previous  approaches  and  rely  on  the  properties  of  the  Interconnection  of 
sub-systems  which  Is  appropriate  In  the  study  of  large  scale  systems. 

2.  We  have  demonstrated  [32]  that  the  spacing  parameter,  which 
appears  in  many  stochastic  approximation  Identification  algorithms  ([K]  and 
[L],  for  example)  Is  unnecessary.  Those  algorithms,  written  as 

jJ(k+£)  »  e(k)  +  ... 

where  A  Is  the  spacing  parameter,  and  k  Is  chosen  as  an  Integer  multiple 
of  A,  should  be  written  as 

Sjk+1)  ■  9{k)  +  ... 

where  k  ■  0,  1,  ...  . 
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IV.  PROFESSIONAL  PERSONNEL  ASSOCIATED  WITH  RESEARCH  EFFORT 


1.  Prof.  Jerry  M.  Mendel,  Principal  Investigator,  1975-1979. 

2.  Prof.  Keith  Glover,  Senior  Investigator,  1975  and  part  of  1976; 
Visiting  Research  Scholar,  Summer  1977. 

3.  Prof.  Leonard  Silverman,  1976. 

4.  Dr.  A.  Ivar  B.  Gustavs son,  Lund  Institute  of  Technology, 

Division  of  Automatic  Control,  Lund,  Sweden;  Visiting  Research 
Scholar,  Summer  1976  (one  month). 

5.  Dr.  H.D.  Washburn,  Ph.D.  student  who  received  Ms  degree  In  June 
1977.  Thesis  title  Is  “Multistage  Estimation  and  State  Space 
Layered  Media  Models." 

6.  Dr.  F.  Hablbl-Ashrafl,  Ph.D.  student  who  received  his  degree  In 
Jan.  1979.  Thesis  title  Is  "Estimation  of  Parameters  In  Loss¬ 
less  Layered  Media  Systems." 

7.  Dr.  F.  Amlnzadeh,  Ph.D.  student  who  received  his  degree  In  June 
1979.  Thesis  title  Is  "Non-Normal  Incidence  State  Space  Models 
for  Layered  Media  Systems." 

8.  D.  Kuan,  J.S.  Lee,  M.  Shiva,  M.  Chan,  and  S.  Kundu,  Ph.D.  students. 

9.  J.  Biddle,  K.  Saedlnla  and  A.  Yu,  computer  programmers. 
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V.  INTERACTIONS 


The  following  papers  which  are  listed  In  Section  III  were  presented 
at  major  technical  conferences:  [2],  [4],  [8],  [9],  [10],  [11],  [12], 
[16],  [17],  [18],  [19],  [20],  [21],  [22],  [24],  [26],  [30]. 

Additionally,  Dr.  Mendel  presented  the  following  talks,  which  were 
based  (all,  or  in  part)  on  research  supported  by  this  grant: 

1.  "Multistage  Least-Squares  Parameter  Estimation:  An  Approach  to 
Modeling  Large  Scale  Systems,"  presented  to  Systems  Engineering/ 
Operations  Research  Seminar  at  the  Unlv.  of  California,  Irvine, 
May  29,  1974. 

2.  "State  Space  Models  of  Layered  Media,"  presented  to  Systems 
Science  Seminar,  Unlv.  of  California,  San  Diego,  Feb.  2,  1977. 

3.  "Estimation  of  Reflection  Coefficients  for  Lossless  Layered 
Systems,"  Unlv.  of  Houston,  E.E.  Dep't.  Seminar,  March  S,  1979; 
EE  Systems  Seminar,  USC,  April  16,  1979;  and  Eindhoven  Unlv.  of 
Technology,  Eindhoven,  The  Netherlands,  Oct.  2,  1979. 
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VI.  NEW  DISCOVERIES  AND  SPECIFIC  APPLICATIONS  STEMMING  FROM 

THE  RESEARCH  EFFORT 

No  patents  Mere  obtained.  The  folloMlng  represent  new  discoveries: 
(1)  generalization  of  Frledland's  [A]  bias  estimation  technique  to  parti¬ 
tioned  dynamical  systems  [6*  7];  (2)  state  space  models  for  layered  media 
systems  [10,  11,  12,  13,  22] i  (3)  a  state  space  Bremmer  series  decomposi¬ 
tion  [11,  14,  29];  and  (4)  a  maximum-likelihood  procedure  for  estimating 
both  the  reflection  coefficients  and  travel  times  for  a  lossless  layered 
media  system  [23,  28]. 

All  of  our  work  described  In  Section  II. C  Is  applicable  to  reflection 
seismology  for  oil  and  gas  exploration. 
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